Here we will recreate the ZooMSS model (version 2) in Heneghan et al. (2020 in preparation) using mizer.
We begin with some setup of required packages.
#get required packages
library(devtools)
#most up to date master branch of mizer
#install_github("sizespectrum/mizer")
#install_github("astaaudzi/mizer-rewiring", ref = "temp-model-comp")
#documentation here:
#https://sizespectrum.org/mizer/dev/index.html
library(mizer)
require(tidyverse)
#remotes::install_github("sizespectrum/mizerExperimental")
library(mizerExperimental)
Next, let’s load in and check the data used in the ZooMSS model. TODO: add in some biomass data
#read in group data
groups <-read.csv("data/TestGroups_mizer.csv")
groups$w_min <- 10^groups$w_min
groups$w_inf <- 10^groups$w_inf
groups$w_mat <- 10^groups$w_mat
#groups<-merge(groups,dat,by.x="species",by.y="group",all=T)
# have a look at plot
# plot(groups$w_inf,groups$biomass.tperkm2,xlab="log Maximum Weight [g]", ylab=" log Total Biomass", log="xy",col="blue",pch=16)
# text(groups$w_inf,groups$biomass.tperkm2,labels=groups$species,cex=0.5)
# could plot the paramter allometries here to explore
Next let’s read in the parameters from ZooMSS.
#read groups again for southern ocean model, this time subsetting key groups
groups <-read.csv("data/TestGroups_mizer.csv")
groups$w_min <- 10^groups$w_min #convert from log10 values
groups$w_inf <- 10^groups$w_inf
groups$w_mat <- 10^groups$w_mat
groups$h <- 10000000000 # should be Inf, but breaks. Massive value still works out to effectively unlimited feeding as allowed in ZooMSS
#groups <- readRDS("data/groups.RDS")[-1,]
#check fails these tests:
#groups$w_mat25 >= groups$w_mat #note we don't have w_mat25. Not sure if it's necessary/desirable. Probably not.
groups$w_mat <= groups$w_min
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
groups$w_inf <= groups$w_mat
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
# read interaction matrix
# get the interaction matrix - actually I think we can leave this out. Default is all 1s, which is the same as in ZooMSS. Included for completeness; it may be useful in future to keep this in.
theta <- readRDS("data/zoomss_inter.rds")[,-1]
#[-1,-1]
We will pass these parameters to mizer to set up a new multispecies model.
TODO: adjust parameters here.
ID <- 1 #index of environmental data to choose
enviro_row <- readRDS("data/envirofull_20200317.RDS")[ID,]
#set up the fixed phyoplankton spectrum
phyto_fixed <- function(params, n, n_pp, n_other, rates, dt = 0.01, phyto_int = enviro_row$phyto_int, phyto_slope = enviro_row$phyto_slope,...) {
n_pp <- 10^(phyto_int)*(params@w_full^(phyto_slope)) #returns the fixed spectrum at every time step
n_pp[params@w_full>params@resource_params$w_pp_cutoff] <- 0
return(n_pp)
}
mf.params <- newMultispeciesParams(species_params=groups,
interaction=NULL, #NULL sets all to 1
no_w = 178, #number of zoo+fish size classes;
kappa = 1e15, #set huge to recreate unbounded feeding on phyto in ZooMSS
min_w_pp = 10^(-14.4), #minimum phyto size. Note use -14.4, not -14.5, otherwise it makes an extra size class
w_pp_cutoff = 10^(enviro_row$phyto_max), #maximum phyto size
r_pp = 20, #Coefficient of the intrinsic resource birth rate, set large for now
n = 0.7, #The allometric growth exponent used in ZooMSS
z0pre = 1, #external mortality (senescence)
z0exp = 0.3,
resource_dynamics = "phyto_fixed",
RDD = constantRDD(species_params = groups) #first go at this
#pred_kernel = ... #probably easiest to just import this/pre-calculate it, once dimensions are worked out
)
## Note: No ks column so calculating from critical feeding level.
## Note: Using z0 = z0pre * w_inf ^ z0exp for missing z0 values.
#checks
length(which(mf.params@initial_n_pp>0)) == length(seq(-14.5,enviro_row$phyto_max, by = 0.1))
## [1] TRUE
Now do some fiddling to make the new MizerParams object match the ZooMSS parameters. Note: this chunk is adapted from fZooMSS_setup.R, found at https://github.com/MathMarEcol/ZooMSS/.
setZooMizerConstants <- function(params, Groups, input_params){
#### CALCULATES CONSTANT BITS OF THE MODEL FUNCTIONS FOR EACH GROUP
SearchVol <- getSearchVolume(params)
M_sb <- getExtMort(mf.params)
ZSpre <- 1 # senescence mortality prefactor
ZSexp <- 0.3 # senescence mortality exponent
pred_kernel <- getPredKernel(params)
prey_weight_matrix <- matrix(params@w_full, nrow = length(params@w), ncol = length(params@w_full), byrow = TRUE)
pred_weight_matrix <- matrix(params@w, nrow = length(params@w), ncol = length(params@w_full))
for(i in 1:nrow(params@species_params)){
## Senescence mortality
if(params@species_params$Type[i] == "Zooplankton"){
M_sb[i,] <- ZSpre*(params@w/(10^(params@species_params$w_mat[i])))^ZSexp
M_sb[i, 10^(params@species_params$w_max[i]) < params@w] <- 0
M_sb[i, 10^(params@species_params$w_mat[i]) > params@w] <- 0
}
if(params@species_params$Type[i] == "Fish"){
M_sb[i,] <- 0.1*ZSpre*(params@w/(10^(params@species_params$w_mat[i])))^ZSexp
M_sb[i, 10^(params@species_params$w_max[i]) < params@w] <- 0
M_sb[i, 10^(params@species_params$w_mat[i]) > params@w] <- 0
}
### Search volume
SearchVol[i,] <- (params@species_params$gamma[i])*(params@w^(params@species_params$gamma[i]))
SearchVol[i, 10^(params@species_params$w_max[i]) < params@w] <- 0
SearchVol[i, 10^(params@species_params$w_min[i]) > params@w] <- 0
### Predation Kernels
if(is.na(params@species_params$PPMRscale[i]) == FALSE){ # If group has an m-value (zooplankton)
# Calculate PPMR for zooplankton, which changes according to body-size (Wirtz, 2012)
D.z <- 2*(3*params@w*1e12/(4*pi))^(1/3) # convert body mass g to ESD (um)
betas <- (exp(0.02*log(D.z)^2 - params@species_params$PPMRscale[i] + 1.832))^3 # Wirtz's equation
beta_mat <- matrix(betas, nrow = length(params@w), ncol = length(mf.params@w_full))
# Calculate feeding kernels
pred_kernel[i, , ] <- exp(-0.5*(log((beta_mat*prey_weight_matrix)/
pred_weight_matrix)/params@species_params$FeedWidth[i])^2)/
sqrt(2*pi*params@species_params$FeedWidth[i]^2)
# The feeding kernel of filter feeders is not expected to change much with increasing size so we fix it here
# if (param$fixed_filterPPMR == TRUE){
if(i == 3){
pred_kernel[i, , ] <- matrix(pred_kernel[i,44,], nrow = length(params@w), ncol = length(mf.params@w_full), byrow = TRUE)
}
if(i == 8){
pred_kernel[i, , ] <- matrix(pred_kernel[i,61,], nrow = length(params@w), ncol = length(mf.params@w_full), byrow = TRUE)
}
# }
} else { # If group does not have an m-value (fish)
beta_mat <- matrix(params@species_params$PPMR[i], nrow = length(params@w), ncol = length(params@w_full))
# Calculate feeding kernels
pred_kernel[i, , ] <- exp(-0.5*(log((beta_mat*prey_weight_matrix)/
pred_weight_matrix)/params@species_params$FeedWidth[i])^2)/
sqrt(2*pi*params@species_params$FeedWidth[i]^2)
}
}
#temperature effect
temp_eff <- matrix(2.^((input_params$sst - 30)/10), nrow = length(params@species_params$species), ncol = length(params@w))
M_sb <- temp_eff * M_sb # Incorporate temp effect on senscence mortality
params <- setExtMort(params, z0 = M_sb)
params <- setSearchVolume(params, SearchVol)
params <- setPredKernel(params, pred_kernel)
}
setZooMizerConstants(mf.params, groups, enviro_row)
Try running it:
sim <- project(mf.params, t_max=5,dt = 0.01)
plot(sim)
plotlyBiomass(sim)
#plotlyGrowthCurves(sim,species="macrozooplankton")
plotlyFeedingLevel(sim)
# feeding level satiation for some groups, except for the seabirds
# macrozooplankton - they are not growing enough,why?
#tuneParams(mf.params)
plotlyGrowthCurves(sim,percentage = T)
plotlySpectra(sim)
—–Old stuff only below this line—–